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In this letter, for the first time, fuU three-dimensional FDTD method to solve the 
Dirac equation is described. The Finite DiflFerence-Time Domain (FDTD) method is a 
fast growing numerical method originally introduced by Kane Yee [T] to solve Maxwell's 
equations. In the last two decades the method was intensively developed primarily in the 
field of electrodynamics [2 [3l [H [5] , but was, at a much lower scale, also extended to other 
fields of applications such as acoustics and elastodynamics. When applied to solving 
Maxwell's equations, the FDTD method is relatively simple and numerically robust, 
has almost no limit in the description of geometrical and dispersive properties of the 
material, and is appropriate for the computer technology of today. As an example, we 
have applied the FDTD method to calculate the exposure of complex biological tissues to 
non-ionizing ultra-wide band (UWB) radiation using high-resolution description of the 
geometry and realistic physical properties of exposed material over a broad frequency 
range [6l |7j . 

In the case of electrodynamics, the FDTD method was able to solve problems with 
complexity that was far beyond allowing analytical solutions and was fundamental to 
the advancement of electrical engineering [5j. In the same sense, we expect that the 
application of FDTD method in quantum mechanics, in this particular case for the 
solution of the Dirac equation, will become a stepping stone for the advancement of 
modern physics. 

Similarities between Maxwell's equations and the Dirac equation are obvious if the 
Dirac equation is written, in a standard notation, as a system of two first-order equations 

m 

h d 

^ cdr 

H^a- V-«^^]$(^) = -mc$(^). (1) 

Two component wave functions and couple in these equations as the electric 
E and magnetic B fields couple in Maxwell's equations [8J. 

While the FDTD scheme can be applied to Eq. Q, it is applied here to the form 
originally written by Dirac. In the case when the electromagnetic field described by 
the four-potential A'^ = {Ao{x)^A{x)} is minimally coupled to the particle, the Dirac 
equation can be written as 

ih— = {Hfree + Hint)^, (2) 



where 



Hfree = —icHa ■ V + (3mc^ , (3) 
Hint = -ea ■ A + eAo, (4) 



and 



'if{x) = 



^2{x) 



(5) 
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Matrices a and /? are expressed using 2x2 Pauli matrices as and the 2x2 unit matrix 
/[SJ. 

The FDTD schematics to solve Eq. ^ foUows Yee's leapfrog algorithm ^j. Space 
and time are discretized using uniform rectangular lattices of size Ax, Ay and Az, and 
uniform time increment At. Any function /(x, t) becomes /(iAx, jAy^ kAz^ nAt) = 
fl^j^k^ where i,j,k, and n are integers. Partial derivatives are expressed numerically with 
second order accuracy as follows: 



(6) 



• space derivative in the x-direction at fixed time nAt 

f) -fn fn fn 

rv^ . 

dx 2Ax 
Analogously for derivatives in y- and z-directions. 

• time derivative at fixed position /(iAx, jAy, kAz) 

~dt 

In the case of electrodynamics the electric field E at the time n — 1/2 is used to calculate 
magnetic field H at the time n, which in turn is used to calculate the electric field E at 
the time n + 1/2, and so on. In the case of the Dirac equation the wave functions and 
at the time n — 1/2 are used to calculate the wave functions and at the time 
n, which are then used to calculate the wave functions and at the time n + 1/2, 
and so on. The same numerical requirements as in the case of electrodynamics were 
followed. Numerical dispersion and stability criteria developed for the electromagnetic 
case [5j were analogously applied to solving the Dirac equation. The Courant stability 
condition was imposed to the time step At 

At < ^ =, (8) 

~ c^(Ax)-2 + (A?/)-2 + (Az)-2' ^ ^ 

and, as a consequence of the Nyquist sampling limit [3j, the rectangular lattices size 

Ax - Ay - Az < A/2. (9) 

A is the wavelength of the plane-wave solution. In practical applications the size of the 
lattice is typically between A/ 10 and A/20. In the case of Zitterbewegung^ the spatial 
oscillation of the wave packet due to an interference between the positive and negative 
energy components [8j , one more condition was imposed to the lattice sizes 

h 

Ax - Ay - Az < . (10) 

^ 2mc ^ ^ 

Boundary conditions appropriate for extending the finite computational domain to 
infinity are applied by extrapolations similar to Liao extrapolation [9j. The physical 
solutions inside the computational domain are extended outside the domain using 
Newton's interpolation polynomials [5]. 

The updating difference equation for was obtained by solving the algebra of Eq. 
([2]), applying Eq. ^ and ([T]), and simplifying using Ax = Ay = Az 
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^1 C-{I,J,K) ^ ^^'^'^J 



+ *^(/ + 1, J, - *^(/ - 1, J, - J+1,K) 

eAt 

- ^^(/, J - 1, K))] + J, (/, J, X) 

- J, K)^2{I. J. K) + J, J, X)], (11) 

where J, X) = 1 + i%[mc^ + eAo(^^ ^)]- Updating equations for ^^2, ^3, 

and were constructed in a similar way. The dynamics of a Dirac electron can be 
now studied in an environment described by any four-potential regardless of its 
complexity and time dependency. As a result of limited space in this letter only few 
cases were chosen. 

Following the analytical solution of K. Huang [TO], we studied the model- 
independent motion of a free Dirac electron acquired from the Dirac equation itself. 
It is important to stress that, as a consequence of the Dirac equation being linear in 
d/dt^ the entire dynamics of the electron is defined by its initial wave function only. 
As in the case of Maxwell's equations, no initial velocity needs to be specified. The 
dynamics of the wave packet studied here was defined by its initial wave function 

/ 1 \ 





^^(x,0) = N\ 



2E 



{pi-ip2)c 



_ x-x I tp-x 

e ^-l ' , (12) 



where is a normalizing constant, A = [(27r)^/^Xo]~'^/^. Eq. (12) represents a wave 
packet whose initial probability distribution is of a normalized Gaussian shape with its 
size defined by the constant Xq. Its spin is pointed along the z-axis and its motion is 
defined by the values of Pi,P2, and J93. In the "single particle interpretation" of the 
Dirac equation, if we choose = = and 7^ the wave packet should move in 
+x-direction. This is not the case. Because of the localization of the wave packet and 
limitation of the direction of the spin, the initial condition in Eq. ( [T2| ) contains also a 
significant negative energy component moving in the -x-direction [lOj. Dynamics of this 
wave packet, for = 18.75 MeV/c and = 10"-^-^ m, is shown in Fig. [ij 

Another aspect of the motion of the Dirac electron consists of Zitterbewegung ^ 
rapid oscillatory motion around a uniform rectilinear motion attributed to interference 
between positive and negative energy states f8]. The FDTD, being a time-domain 
method, is an excellent tool to study the properties of Zitterbewegung. Few results 
are presented here. Fig. [2] shows the time dependence of the position of the center 
of probability in the x-y plane of the fraction of our wave-packet moving in the -x- 
direction at the beginning of motion. The amplitude of the motion in this particular 
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Figure 1. Top view of the four different stages of motion of the wave packet initiahzed 
by Eq. (12) in the x-y plane. The initial wave packet splits in two, moving in both 
the positive and negative x-directions, corresponding to positive and negative energy 
solutions. It is important to notice details of electron dynamics through changing of 
the shape of the wave packets into elliptical as a result of relativistic suppression of 
the wave packet spreading in the direction of motion for a massive particle pT]. The 
full animations can be accessed on-line [12.. 
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Figure 2. Time dependence of the position of the center of probabihty of the fraction 



of the wave-packet defined by Eq. ( 12 ) moving in the -x-direction at the beginning of 
motion. 



case is 1.45 x 10~^^ cm, of the order of h/mc as predicted [8j. Fig. [3] shows the x-position 
of the fraction of the wave-packet moving in the +x-direction. Zitterbewegung exists at 
the beginning of the motion, at the time when positive- and negative-energy components 
of the wave packets shown in Fig. [T] overlap. It is absent as the wave packets separate. 
Lack of space in this letter prevents discussion of Zitterbewegung-ieldited violent and 
rapid oscillations of the velocity of the electron. 

The frequency of Zitterbewegung can easily be obtained. Fitting the probability 
function, as shown in Fig. [4| we found the angular frequency of cu = (1.582 ± 0.009) x 
10^1 Hz for pi = 18.75 MeV/c and uj = (1.525 ± 0.009) x lO^i Hz for pi = 187.5 eV/c. 
It is in agreement with expected angular frequency of ^ 2mc^/h and shows little 
momentum dependence. 

Finally, we discuss briefly the well-localized state described, for example, by 



*(f,0) = N 



_ x-x 
4x'i, 






Snapshots of the dynamics of the probability densities of this state for xq 



(13) 



10 



-14 



simulating localization in a volume of the size of the atomic nucleus, are shown in Fig. [5) 
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Figure 3. X-position of the fraction of the wave-packet shown in Fig. [T] moving in 
the +x-direction. The Zitterbewegung which exists at the beginning of the motion 
(top) is absent after positive- and negative-energy components separate (bottom). 
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Figure 4. Fit of the probabihty function used to determine the frequency of 
Zitterbewegung. 

FDTD computation shows that the angular distribution of probabihty density is 
determined by spherical harmonic yoo, the density |^^2p is 0, the angular distribution of 
|^^3p is determined by spherical harmonic l^io, and |^^4p by The radial dependence 
of l^^ip is determined by the spherical Bessel functions jo, and |^^3p and |^^4p by 
ji. The time dependent behavior of the total probabilities is shown in Fig. |6| After 
the first shock, attributed to Zitterbewegung resulting from the initial condition, total 
probabilities acquire the values of Probl : Prob2 : ProbS : ProbA = 1/2:0:1/6:1/3. 
The ratios of those probabilities result from the Clebsch-Gordon coefficients which are 
part of the relation between spherical harmonic spinors ^jim and spherical harmonics 
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Figure 5. Three different stages of motion of the wave packet in the x-z plane 
initiahzed by Eq. (13). Top-row is probabihty density mid-row |^3p, and 



bottom-row |^4| . The associated animations can be accessed on-hne [12]. 



Yim p^fT^ . (Those ratios persist for more localized states, for example for xq = 10~^^ m, 
but, as expected, change for less localized states.) The angular and radial distribution 
of the probability densities and the values of the total probability determined by the 
FDTD method show that the state with initially localized position and spin, corrects 
its state in a very short time (Fig. [6]) and behaves as a spherical wave. 

To conclude, for the first time full three-dimensional Finite Difference Time Domain 
(FDTD) method was developed to solve the Dirac equation. In this paper the method 
was applied to the dynamics of a free Dirac electron, comparing some of the results to 
the analytical solutions. Forthcoming work will study the dynamics of the Dirac electron 
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Figure 6. Time dependence of the total probabilities of the wave packet defined by 
Eq. (p). 



in potential fields. The impact of FDTD method in electrodynamics warranted that the 
first chapter of the book by A. Taflove and S. C. Hagnes [5j be titled "Electrodynamics 
Entering the 21st Century" . We hope that the FDTD method will have the same impact 
on better understanding, advancing, and applying modern physics. 
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